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INTRODUCTION 



o 

u 

\ A self-gravitating system (SGS) is a system where many particles interact via the gravi- 

tational force. When we shall explain a distribution of SGS in phase space, the Boltzmann- 
Gibbs statistical mechanics is not useful. This is because the statistical mechanics is con- 



structed under the condition of the additivity of the energy: as well known, the total energy 
of several SGSs is not equal to the sum of the energy of each system. In fact, SGS does not 

> 

CN have a tendency to become a state characterized with a temperature. 

m 

■ If we use the statistical mechanics assuming that the state of the SGS with an equal 

QQ . mass m becomes isothermal with the temperature T and that the particles of the system 

o ■ 

are distributed spherically symmetrically, what kind of distribution can be obtained? Then, 
the structure in phase space can be determined using the Maxwell-Boltzmann distribution. 



^ \ For example, the number density at a radial distance r in real space is given by 

nMB{r) oc e '^b^ ^ > ^ (1) 

where $(r) is the mean gravitational potential per mass generated by this whole system at 
r and k-Q is the Boltzmann constant. This potential should satisfy a relation with number 
density by the Poisson equation, 

— rV + ^ = 47rGmnMB(r) , 2 

dr^ r dr 
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where G is the gravitational constant. A special solution to Eq. ([T]) and this Poisson equation 
is ny[^{r) = k-^T , known as the singular isothermal sphere However, this 
solution has two problems: infinite density at r = and infinite total mass. Even though 
we solve the equations with a finite density at r = 0, the solutions behave oc T at a large 
r, and so we cannot avoid the infinite total mass problem. In either case, the solutions are 
unrealistic. 

Of course, real examples of SGS in the universe, e.g., globular clusters and galaxies, have 
various structures with a finite radius. As for most globular clusters, it is known that their 
number densities in real space have a fiat core and behave as a power law outside this core. 
King interpreted these profiles by introducing the new distribution function in phase space, 
known as the lowered Maxwellian; 

( g-/3(i?-™<i.o _ 1 for E < 
f{r,v)^{ - (3) 

I for E > m^t , 

in which E is the total energy of a particle belonging to a globular cluster. This distribution 
becomes zero when the total energy is greater than m^t, and so $t can be understood as a 
potential energy per mass at the surface of the globular cluster. Because the velocity of the 



particle must be in the range of < v < a/2 — $(r)}, the number density rixuiT") can 
be obtained integrating f{r,v) as 



•i/2{$t-$(r)} 

riKMir) (X I dvATTV^f{r,v) . (4) 



Moreover, using a dimensionless potential W{r) = — m/3 {<l>t — <l>(r)} and integrating by 
parts, the number density becomes 



nKuir) oc e^^'') 



/ dCe-^C^' . (5) 

^0 



As mentioned before, the potential energy and the number density has a relation through 
the Poisson equation. Thus, W{r) must satisfy the following equation: 

dW(r) ^ 2 dW{r) _ 9 UKuir) 

where a = a/9/ {47rG'm^/3nKM (0) } corresponds to the core radius. The number density 
satisfying Eqs. ([5]) and (|6]) can be calculated numerically as shown in Fig. [1] This is called 



the King model When W{0) is larger than about 5, the number density around the 
origin can be represented by the following approximation: 

1 



riKuir) oc 

which is shown as the red curve in Fig. [H 



(l + rVa2)3/2 



0.01 



0.0001 



1.x10"^ 

0.01 



(7) 
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FIG. 1: Black curves denote the number density of the King model raKM(^) as a function of the 
radius normalized with the core radius, r/a, for several VF(0). As the curve changes from left to 
right, VF(0) gets larger from 2 to 10 in steps of 2. The red curve means the approximated formula, 
1/(1 + r2/a2)3/2, for larger W{0). 



Since King put forward this model, this number density has been applied to fitting for 
the surface brightness of many globular clusters, for example as in Ref. |3|-[7|, that is, most 
exponents of power law outside the core of globular clusters are similar to —3 which cannot 
be explained by the model with the isothermal assumption. But, it is not easy to see what 
kind of dynamics occurred in the system, because his procedure was done to the distribution 
function in the steady state. 

So, we will construct a theory which can explain the dynamics toward such a special steady 
state described by the King model especially around the origin. The idea is to represent an 
interaction by which a particle of the system is affected from the others by a special random 
force described by a position-dependent intensity noise, in other words the multiplicative 
noise, that originates from a fluctuation only in SGS. That is, we will use a special Langevin 
equation, just as the normal Langevin equation with a constant-intensity noise, in other 
words the additive noise, can unveil the dynamics toward the steady state described by the 
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Maxwell-Boltzmann distribution. However, we cannot introduce the randomness into the 
system without any evidence. Then, we must confirm that each orbit is random indeed. Of 
course, it is impossible to understand orbits of stars in globular clusters from observations. 
Thus, we must use numerical simulations. 

From the numerical simulations of SGS, the ground that we can use the random noise 
becomes clear. The special Langevin equation includes additive and multiplicative noises. 
By using this stochastic process, we derive that non-Maxwell-Boltzmann distribution of SGS 
especially around the origin. The number density can be obtained through the steady state 
solution of the Fokker-Planck equation corresponding to the stochastic process. We exhibit 
that the number density becomes equal to the density profiles around the origin, Eq. (^^, 
by adjusting the friction coefficient and the intensity of the multiplicative noise. 

Moreover, we also show that our model can be applied in the system which has a heavier 
particle (5-10 times as heavy as the surrounding particle). The effect of the heavier particle 
in SGS, corresponding to a black hole in a globular cluster, has been studied for long time. If 
the black hole is much heavier than other stars, a cusp of the density distribution appears at 

- lo||. The observations which suggest that intermediate mass black 



the center of a cluster 



hole (IMBH, ~ 10^ — IO^Mq) is in the globular cluster in recent years are accomplished one 



after another 
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17] . Although these studies are very interesting, our model does not treat 



these situations: in our model, the heavier particle is too lighter than IMBH. Our model 
corresponds small globular cluster (10^ stars) with only a stellar black hole (~ 1 — IOMq). 



Here, note that we have reported similar results in our previous letter |ll|. In this paper, 
however, we demonstrate how we executed our numerical simulations. Moreover, a treatment 
for stochastic differential equations becomes precise, and so the analytical result derived by 
a different method changes a little. 

This paper is organized as follows. In Sec. Ill At and Sec. Ill Bt we provide brief explana- 
tions about a machine and a method we used when we did numerical simulations, respec- 
tively. In Sec. Ill CI we investigate number densities derived from our numerical simulations 
where all particles of SGS with a mass m and a particle with a mass M interact via the 
gravitational force. Then, we show the densities are like that of the King model and both 
the exponent and the core radius are dependent on M. In Sec. IIII Al forces influencing each 
particle of SGS are modeled. Then, using these forces, Langevin equations are constructed 
in Sec. IIIIBI Section llH CI makes it clear that the steady state solution of the corresponding 
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Fokker-Planck equation gives the same result with the King modeL In Sec. llVt we discuss 
our results and make the relation between King's procedure and our idea clear. Section |V] 
gives a summary of this work. 

II. NUMERICAL SIMULATIONS OF SGS USING GRAPE 

A. GRAPE 

SGSs require quite long time for relaxation. Furthermore, because only attractive force 
is exerted on particles in SGS and the gravitational potential is asymptotically flat, we 
must compute interaction of all particle pairs. When we treat N particles, the computa- 
tion of interaction becomes 0{N'^) by direct approach. By these reasons, we require huge 
computation for numerical simulation of the evolution of SGS. 

For time evolution of SGS, many improvements of algorithm and hardware have been 
carried out. First, we consider integrator for simulation. For long-time evolution, both 
the local truncation error and the global truncation error are noticed. These error occur 
deviation of the conservation physical quantities such as total energy. For compression 
of the global truncation error, symplectic integrator has been developed. The symplectic 
integrator conserves the total energy for long-time evolution. We apply 6th-order symplectic 
integrator for the time evolution of SGS. Secondly, we apply special-purpose processor for 
the computation of the interaction. Most of the computation of the time evolution in SGS is 



2-body interaction. As special-purpose processors, GRAPE system has been developed jl8 |. 
GRAPE system can compute 2-body interaction from position and mass of particles quickly. 
In our study, we apply GRAPE-7 chip, which ^consists of Field-Programmable Gate Array 
(FPGA) for computation of the interaction [19]. GRAPE-7 chip implements GRAPE-5 
compatible pipelines jsSj. The performance of GRAPE-7 chip is approximately 100 GFLOPS 
and is almost equal to a processor of present supercomputers, but the energy consumption 
of the chip is only 3 Watts. Using sophisticated integrator and special-purpose processor, 
we have analyzed time evolution of SGS. 
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B. Symplectic integrator 

For time evolution, we must choose reasonable integrator for simulation. For long-time 
evolution, not only the local truncation error but also the global truncation error is noticed. 
For example, 4th-order Runge-Kutta method has been applied for time evolution of physical 
systems l^l]. Although its local truncation error is O ((At)^), because its error accumulates, 
the global truncation error increases during time evolution. For example, we apply 4th-order 
Runge-Kutta method for the harmonic oscillation. 

Using exact solutions, the orbit of the harmonic oscillation in the phase space draws a circle. 
Of course, the total energy is conserved. On the other hand, when we apply 4th-order 
Runge-Kutta method for time evolution, the total energy is decreased monotonically. 



Hit) = - 



l(Atf + 0((Atf) 



[P' + <f) ■ (9) 



The orbit of the harmonic oscillation in the phase space draws a spiral and it converges to 
origin = ^7 = 0). When we consider long-time evolution, 4th-order Runge-Kutta method 
is not reasonable integrator. If Hamiltonian of physical system is given, we can apply 



symplectic integrator which based on canonical transformation [22|-[2J, |26|| . This method 
suppresses increase of the global truncation error. In generic case, the symplectic integrator 
is implicit method. If Hamiltonian is divided to coordinate parts and momentum parts, 
the integrator becomes explicit method. The procedure of low-order integration becomes 
easy more than Runge-Kutta method. The simplest integrator is called "leap- flog method" 
(2nd-order integrator). 

+ — j = p{t) + —p{x{t)) , (10) 
x{t + At) = x{t) + At ■ p (^t + , (11) 
p{t + At) = p (^t + + {x{t + At)) . (12) 
Using leap-flog method for the harmonic oscillator, the following equation is satisfied. 

+ g^) + ^pg = const. . (13) 
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Therefore the orbit in the phase space draws an oval. The deviation from the exact solution is 
suppressed. To suppress the local truncation error, higher-order symplectic integrators have 
been developed. We apply 6th-order symplectic integrator for time evolution of SGS [25 1. 

Pi = + c,Atp{q,^i) (1 < i < 8) , (14) 
qj = pj^i + djAtp, (1 < J < 7) , (15) 

where po = p(t), go = q{t),ps = pit + At),qj = q(t + At). The coefficients Cj, and dj are 
shown in Tab. [B 



i 






1 


0.39225680523878 


0.78451361047756 


2 


0.510043411918458 


0.0235573213359357 


3 


-0.471053385409757 


-1.17767998417887 


4 


0.06875316825252 


1.31518632068391 


5 


0.06875316825252 


-1.17767998417887 


6 


-0.471053385409757 


0.0235573213359357 


7 


0.510043411918458 


0.78451361047756 


8 


0.39225680523878 





TABLE I: Coefficients of 6th-order symplectic integrator (Solution A in 



The symplectic integrator conserves the total energy and the symplectic structure in 
generic cases. When we use n-th order symplectic integrator, the local truncation error 
of the total energy becomes 0((At)"+^). Furthermore, the global truncation error is not 
accumulated 

In SGS, because the interaction at zero distance diverges, the local truncation error would 
diverge in long-time evolution. For avoidance of this divergence, some kind of softening 
parameter has been introduced to gravitational interaction. When the nature of the pure 



gravity is analyzed, the regularization procedure of interaction is required 



28|, 



29|. 
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C. Steady number density in numerical simulation 

Now, we investigate the steady number density (SND) of the SGS with a mass m including 
a particle with a mass M by numerical simulation. In particular, we show that SNDs have 
a core and behave as a power law outside the core. 

The system is composed of = 10000 particles. At t = 0, all velocities of the particles 
are zero and they are distributed by no(r) oc (1 + r^/op^)^^/^ (0 < r < 4ap), which is the 
density in real space of Plummer's solution In this SGS, we put another particle with 



a mass M in the origin at t = 0. We shall change the mass as M/m = 1, 5, and 10. 
Throughout this paper, we adopt a unit system where the core radius of Plummer's solution 
Op, the initial free fall time tfj, and the total mass ■ m are unity. 

We started the numerical simulation under these conditions. For dynamical evolution, we 
used GRAPE-7 at Center for Computational Astrophysics, CfCA, of National Astronomical 
Observatory of Japan. For the computation of gravitational force, we applied Plummer's 
softening: the potential energy between the ith and jth particles separated by a distance 



Vij is —Gm'^/wrij^+e^, where Eg is the softening parameter. We set Eg = 10 ^. For time 



evolution, we used 6th-order symplectic integrator 25|. The time step for the simulation is 
defined as At = 10~^. We carried out simulations until t = 100 tfj. During simulations, the 
error in total energy was less than 0.1%. 

First, most particles collapse into the origin within several tff. After approximately 
20 tff, the distribution becomes stable and the system reaches the steady state. Of course, 
we can confirm whether the system becomes steady or not from the profile of the number 
density. However, furthermore we also focus on the number of particles inside a sphere. 
Figure [2] shows the change of the number inside the sphere with a radius 1 in time. During 
the collapse, the number becomes large. After that, the number decreases, which means 
that many particles with positive energy evaporated from inside of the sphere, and so the 
number becomes about 6000 on average. For other radii, similar changes of the number in 
time can be seen. 

SND is calculated by taking the time average during the steady state. In Fig. |3l we show 
the logarithm of SND as a function of logr for M/m = 1,5, and 10. For each M, the SND 
has a core and behaves as a power law at r larger than the core radius. Here, we fit SNDs 



around the core by nfit{r) = C/(l + r"^ /a^Y. The results are summarized in Tab. [TTl For 
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FIG. 2: Change of the number of particles inside a sphere with a radius 1 in time for M/m = 1, 
5, and 10. 

M/m = 1 and 5, k ^ 3/2, which is similar to the exponent of the King model. The density 
at the origin C increases as M increases, which is simply understood to be a result of many 
particles being attracted by the heavier particle. 



M/m a X 10^ 



C X IQ-^ 



6.20 ±0.22 1.46 ±0.03 7.06 ±0.14 



6.06 ±0.18 1.46 ±0.02 7.57 ±0.13 



10 5.68 ±0.14 1.43 ±0.02 8.13 ±0.12 



TABLE II: Best-fit parameters of the function C/(l ± r'^/a^)'^ for steady number densities shown 
in Fig. [3] 



III. SIMPLE MODEL 

A. Forces acting on each particle of SGS 

As shown in the last section, SND is the King-like profile even though the system includes 
the heavier particle. In this section, in order to explain these results and derive this non- 
Maxwell-Boltzmann distribution around the origin, we demonstrate a simple model based 
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FIG. 3: Logarithm of the steady number density derived from our numerical simulation, n(r), as 
a function of logr for (a) M/m = 1, (b) M/m = 5, and (c) M/m = 10. In each figure, the dashed 
curve and symbols denote a fitting curve n^t{r) = C/(l + r"^ /a?)'^ and the result of our numerical 
simulation, respectively. 
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on stochastic process, which is quite different from the King modeL 

The reason why stochastic process appears in the SGS is as follows. After the collapse, 
the density around the origin becomes high. Thus, the particles around the region disturb 
the orbits of other particles repeatedly, so that their movements become random 36|]. As 
the time at which this disturbance occurs, we introduce the local 2-body relaxation time 

^rel 



trelir) 



0.065(T r 



(16) 



G'^n{r)m? \u.{l/es) 

where a{r) is the standard deviation of the velocity at r; we adopted \n{l/es) as the Coulomb 
logarithm. 



Figure m shows the logarithm of trei, which is calculated using the a{r) and n(r) during 
the steady state obtained from our numerical simulation, as a function of logr. As expected, 
trei around the origin is short. Our simulation continues after the collapse at about 80 t//, 
which is sufficiently longer than the trei around the core. As radius increases, however, t^ei 
becomes longer than the rest of our simulation time, which means that no stochastic motion 
occurs at a large r. Therefore, note that our model is valid only in the neighborhood of the 
core. 




FIG. 4: Logarithm of the local 2-body relaxation time tj-ei as a function of logr for M/m = 1, 5, 
and 10. 



When constructing our model, the following points are premised: the model describes the 
stochastic dynamics near the steady state and the mean distribution is spherically symmetric. 
As is well known, the gravitational force at r arising from such a spherically symmetric 
system depends only on the particles existing inside a sphere with a radius r, and this 
attractive force acts along the radial direction. In other words, this mean force —F{r) is 
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the gradient of the mean potentiah — F(r) = —md^{r) /dr (< 0). Indeed, Mmr^Q F{r) = 0. 
Hence, F[r) can be expanded around the origin as 

F{r) = Cor + 0{r^) . (17) 

It will be clarified later that the lowest exponent must be 1 and the coefficient Cq is related 
to the number density at the origin C as 

Co = 47rGm2C/3 . (18) 

For M/m = 1, we can identify another particle together with the other particles. Contrary 
to this, we must consider the effect of the particle in the case M/m ^ 1. Now, we suppose 
that the heavier particle exists at the origin. Then, the attractive force by this particle 
at r is — t^(r) = —GmM/r"^. We can estimate F{r) around this region as F{r) ~ c^r = 
AirGm'^Cr/S, where we used Eq. ([H]). Thus, ^(r)/F(r) ~ SMr"^ /ATrCm. This ratio 
becomes signiffcant when r ^ 10^^, since C ~ 10^ as shown in Tab. [Tll Therefore, if r is 
smaller than the radius, particles are inffuenced by not only F{r) but also <^{r), so that the 
core disappears. In fact, we have performed a numerical simulation with the heavier particle 
ffxed at the origin, where this result is conffrmed. On the other hand, Miocchi improved the 
King model in order to describe the steady state of a globular cluster including an IMBH and 



reported that the density becomes cuspy as the mass of the black hole increases [3l|]. The 
nature of a globular cluster when a massive black hole is much heavier than the surrounding 
star, have been studied as mentioned in Introduction. In this case, the massive black hole 
stays at the center mostly. Then, a cusp of the density distribution at the center appears. 
Because the heavier particles in our numerical simulation are not very heavy, the particles 
are not trapped at the origin. Therefore, we do not consider the effect of heavier particles 
explicitly and we suppose that the particles inffuence SGS through the density at the origin 
C: as M becomes larger, it attracts more particles and C increases, as shown in Tab. [Tll 
Thus, Cq is an increasing function of M. 

It is natural to consider that the distribution ffuctuates around the mean because of 
the many disturbances. In fact, as shown in Fig. [21 the number of particles existing inside 
the sphere with a radius 1 ffuctuates around the mean value. The ffuctuating part of the 
distribution should not be spherically symmetric, so that this produces forces along not only 
the radial direction, but also other directions. We assume that they are random forces and set 
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their intensity at r 2H{r). In addition to such random forces resulting from the fluctuating 
distribution, a particle at r is expected to be influenced by random forces generated from 
neighboring particles. We set the intensity 2D, which is independent of position. 



B. Langevin equations 

Stochastic dynamics under the above assumptions is described by the following Langevin 
equations in spherical coordinates: the radial direction 



mar + mjr = -F(r) + ^/2H{r)rir{t) + V2D^rit) , (19) 
the elevation direction 



mag + m-fr9 = ^2H{r)7]e{t) + V2Die{t) , (20) 

and the azimuth direction 

ma^ + m7r sin 6^ = ^2H{r)r]^{t) + V2D^^{t) , (21) 

where a^, ag, and are accelerations along those directions; 7 is the coefficient of dynamical 
friction in the low velocity limit, independent of velocity iL In the Chandrasekhar dynamical 
friction formula, the coefficient is more complicated l|, |32| . However, we use the coefficient 
in such a limit, because the density around the core is so high that particles around there 
move slowly. 

Now, we focus on the overdamped limit of these equations, because we have interests in 
the stochastic dynamics near the steady state. In the case of the normal Langevin equation 
with a constant-intensity noise, we only neglect the inertial term. But, as for specm/ Langevin 
equations with noises whose intensity^ depends on a position, the new force —'VH{r)/2m'y 
should be considered additionally [37|]. Thus, the Langevin equations in the overdamped 
limit becomes as follows: 



radial direction : m-fr = -F(r) + J2H{r)r]r{t) + V2D^r{t) —H'(r) , (22) 

2m7 

elevation direction : mjrO = ^/2H{r)rig{t) + \/2DC,g{t) , (23) 



azimuth direction : m7rsin^^0= y'2H{r)ri^{t) + V2D^^{t) , (24) 
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where the prime indicates a derivative with respect to r. The noises in each Langevin 
equation, ^i{t) and rij{t) = r, ^,and (p), are zero-mean white Gaussian and correlate 



only with themselves. Indeed, the correlation function is the Dirac delta function [38|. 

Here, revisit the position-dependent intensity noise. We have introduced such a noise in 
order to represent a random force originating from the fluctuation of the distribution around 
the mean value which yields the mean force —F{r). Thus, the first and the second terms on 
the right-hand side of Eq. fl22|) must denote that the mean force acting along radial direction 
is fluctuating. As a minimal formulation describing this situation, we propose the following 
one: 

-F{r){l-V2'er]r{t)} , (25) 
in which e is a positive constant. This can be realized by setting 

H{r) = eF{rf . (26) 

Note that this fluctuating mean force is the essential feature for SGS. Since the gravitational 
force is a long-range one, each particle is influenced from the whole system. The mean force 
is produced by the mean potential which is decided by the number density in the steady 
state through the Poisson equation. Obviously, this number density determines only the 
mean positions of the particles, and they do not remain stationary at those positions: they 
fluctuate. Then, the mean force also fluctuates, e indicates the extent of fluctuations. If e is 
0, that is, means the mean force does not fluctuate, the stochastic dynamics of each particle 
is governed only by the constant-intensity random force originating from the neighboring 
particles. Then, the Maxwell-Boltzmann distribution is obtained as the steady solution, by 
which the number density of globular clusters cannot be explained as written in Introduction. 

C. Fokker-Planck equation and the asymptotic steady solution around the origin 

From the Langevin equations (!22l) . ( |23l) . and (|24l) . we obtain the Fokker-Planck equation 
governing the spherically symmetric probability distribution function (PDF) P{r,t) 

<9 . D r 92 2 9 1,, lid 



+ i [—Firf + -—F{rf\ P{r,t) , 
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in which we have replaced H{r) by F{r) using Eq. (126|) . Then, the PDF with the Jacobian 
p(r, t) = 47rr^P(r, t) satisfies the following Fokker-Planck equation. 

d , , D ( 8 2]^^ 1 5 ^ . ^ 

+ 7Z^II^-|:-VwV(r,t) (27) 



(7717) 2 1^ dr'^ dr r 

This equation is useful when integrating with respect to r. 

The steady state solution Pst{r) satisfies Eq. fl27j) with the left-hand side zero. By inte 
grating the equation with respect to r, we have 



I 2F{r)F'{r) - -F{rf \ 



Pst{r) = const. . (28) 



(m7)^ r (7717)^ 1^ r j 7717 

Now, we impose the binary condition that Pst(^) = Pst('")/(4vrr^) and the derivative do not 
diverge at the origin. Then, when r — )■ 0, 

Pst(r) = 0{r') , (29) 

and 

lim = lim4vr(2rP,t(r) + r2p4(r)) = , (30) 

r— >-U r— >0 

by which the constant on the right-hand side of Eq. fl28|) is decided and we obtain 

r{D + eF{rf} p'^^{r) = - [rP(r) {2eP'(r) + m-f} - 2 {D + eP(r)2}] pst(r) . (31) 

Thus, if F{r) is obtained, Pst{i^) can also be obtained. Here, F{r) relates with SND, 
?7,st(r), through the following relation, since —F{r) = — m$'(r). 

F'(r) + -F{r) = AnGm'^n.tir) (32) 
r 

Incidentally, the SND can be obtained by multiplying PDF in the steady state by total 
number N: 

nst(r) = iVP,,(r) = . (33) 

Therefore, equation (l32l) can be represented as 

p,. .^2 GNm^p,,{r) 

F (r) + -P(r) = . (34) 
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In short, Psti'f') and F{r) are closely connected with each other through Eqs. (13T|) and 
f p^ . Here, we focus on the asymptotic behaviors of them around the origin, since our model 
is valid around there as mentioned before. Furthermore, due to this approach, we can treat 
them analytically. 

Assume that F{r) can be expanded around the origin with the lowest exponent k as 
follows. 

oo 

F{r) =r''J2 (35) 

1=0 

Substituting this expression into Eq. f l34|) . we find that Pst('") can also be expanded like 

Pst(r) = ^^5^Q(fc + / + 2)r'. (36) 

1=0 



After substituting both Eqs. fl35|) and fl36|) into Eq. (l3T|l . we can obtain 
D + er'' Yl 7^ + / + 1)(A: + / + 2)r' 



r 

X 



1=0 

oo f oo ^ I / oo 

^ Qr' <^ 2er'=-^ ^ c,(fc + s)r^ + 7717 I - 2 i D + er^'^ I ^ Qr' 

/=0 t s=0 J [ \l=0 

k+1 °° 

J2ci{k + l + 2y. (37) 



GNm^ 

1=0 



Firstly, we compare the lowest order terms on the both hand sides of Eq. ( 137|) . so that 
the following relation can be seen: 

Therefore, we can conclude that k = 1. Secondly, compare the next lowest order terms 
proportional to and we get 

and so Ci = 0. Lastly, selecting only terms proportional to from Eq. fl37p . we can find 

from which the following relation can be obtained: 

5 C2 eco^ / m'j 



3 Co -D V 2eco 



l + • (41) 
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Without going into detail, we can see that C3 = by comparison with terms containing r^. 
So, Pst{r) becomes as follows: 

Pst(r) = ^^(3co + 5c2r^)+0(r'') 

3cor^ fl-!^fl + ^V4+0(r^)- (42) 



GNm^ 1 D \ 2eco 



Here, if we set 



2 ^ '^1 
a = — ^ and k = 1 H , (43) 

eco^ 2eco 
Pst('") can be expressed around the origin like 



Pst(r) ~ ^^^2(i + ^2/^2^K ' (44) 



which yields 



Thus, we can derive the number density around the origin of SGS from the model using 
stochastic dynamics. 

The relation flTSl) is easily obtained by setting r = on Eq. fHSl) . that is, 

IV. DISCUSSION 

In this section, we investigate the results derived in the preceding section and understand 
the roles of two noises and the heavier particle in Eq. f HS]) . Additionally, we discuss the 
difference between the King model and our model. 

As in Eq. ( H3|) . the exponent k must be larger than 1, which does not contradict our 
numerical simulation shown in Tab. [Tll In order for Eq. f l45|) to correspond completely to 
the King model, k = 3/2 or 7 = ec^/m = AnGmeC/S must hold. We can regard this 
relation between the friction coefficient 7 and the intensity of the multiplicative noise e as 



a kind of fluctuation- dissipation relation [3J], which usually plays an important role when a 
stochastic process with a constant-intensity noise goes to the equilibrium state described by 
the Maxwell-Boltzmann distribution. 
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The core radius a is proportional to a square root of the intensity of the additive noise 
D owing to Eq. f l43|) . Then, this intensity spreads the region where the density is almost 
constant. This is recognized as the effect of this noise, which makes a system homogeneous 
and isothermal. The existence of the core at globular clusters shows that such a diffusive 
effect does not disappear even for the system with long-range force. In other words, all 
statistical mechanical features observed in a system with short-range force, that is, normal 
system, does not change drastically in SGS and this effect is still universal. Our model 
makes it clear that the special distribution can be obtained just considering the fluctuation 
of mean force. 

Now, let us examine the role of the mass of the heavier particle, M, in this system by a 
naive discussion. As mentioned previously, cq is an increasing function of M. cq exists in 
the denominators of a and k. Then, both values should be reduced when M is increased if 
other parameters are independent of M. These theoretical expectations are consistent with 
our numerical results shown in Tab. [TTl 

How the special distribution (H5l) changes if we do not consider the fluctuating mean 
force? The steady state solution of Eq. ( 1271) with e = is 



Therefore, our result goes to a singular isothermal sphere, as discussed at the beginning of 
this paper, by which the number density of globular clusters cannot be explained. 

Here, we examine the relationship between the King model and our model. King trans- 
formed the distribution function in the phase space in order to avoid a singular isothermal 
sphere. In our model, we introduce the multiplicative noise into the system influenced by 
the mean force and the additive noise whose PDF becomes Maxwellian in the steady state, 
as shown in Eq. (H7j) . so that the non-Maxwell-Boltzmann distribution ( l45l) is derived. In 
short, although these procedures seem to be different, they may have the same meaning at 
least around the origin. However, we emphasize that the stochastic dynamics around there 
near the steady state becomes clear owing to Eqs. (fT9|) . (l20l) . and (I2T]) . 



2 




(47) 
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V. CONCLUSION 

In conclusion, the non-Maxwell-Boltzmann distribution ( H5|) has been obtained using the 
stochastic dynamics with the fluctuating mean force and the additive white noise. This 
number density can be the same as that of the King model around the origin by controlling 
friction coefficient and the intensity of multiplicative noise. Furthermore, our model can 
describe the SGS having a heavier particle. Of course, these results are consistent with 
our numerical simulation. We can say that such a stochastic dynamics occurs behind the 
background of the King model. In short, the diffusive effect, which is represented by the 
additive noise, is universal even in SGS, and it is particular to SGS that the fluctuation of 
the distribution around the mean value producing the mean force makes influence on each 
particle of this system, which our simple model can describe. 

Finally, note that our result is available only in the neighborhood of the origin. Therefore, 
we must derive the density globally by further extended model and investigate the difference 
between the model and the King model. 



Appendix A. DIMENSIONS OF ^/D/eco^ AND m'y/2ecQ 

From now on, [•] represents a dimension of •. Since the correlation function of the random 
noises ^j(t) and rjj{t) = r, 9, and 0) is the Dirac delta function with argument t, 

[e.(^)] = fe(^)]=time-l/^ (A.l) 

Thus, from the expression f l25p whose dimension is a force, we can see that 

[e] = time . (A.2) 

Furthermore, from \/2D^r{t) whose dimension is also a force, the dimension of D can be 
clear like 

[D] = force^ ■ time = mass^ ■ length^ ■ time~^ . (A. 3) 

Owing to Eqs. f lT7|) or f|T8|) . the dimension of Cq equals a force per length: 

[co] = force • length"^ = mass • time~^ . (A. 4) 
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As well known, the dimension of the damping constant, 7, is an inverse of time: [7] = 
time~^. Thereby, 

y 

and 

7717 

eco 
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mass^ ■ length^ ■ time ^ , , , . 

2 -4 = length , (A.5) 

time • mass^ • time 



mass • time 
time • mass • time" 
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